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Abstract 

The current status of computational methods for unsteady 
aerodynamics and aeroelasticity is reviewed. The key features 
of challenging aeroelastic applications are discussed in terms of 
the flowfield state: low-angle high speed flows and high-angle 
vortex-dominated flows. The critical role played by viscous 
effects in determining aeroelastic stability for conditions of in- 
cipient flow separation is stressed. The need for a variety of 
flow modeling tools, from linear formulations to implementa- 
tions of the Navier-Stokes equations, is emphasized. Estimates 
of computer run times for flutter calculations using several com- 
putational methods are given. Applications of these methods 
for unsteady aerodynamic and transonic flutter calculations for 
airfoils, wings and configurations are summarized. Finally, rec- 
ommendations are made concerning future research directions. 

1. Introduction 

In the past decade there has been much activity in the de- 
velopment of computational methods for the calculation of un- 
steady aerodynamics about airfoils, wings and complete vehicle 
configurations. Two key areas of activity have been transonic 
aeroelasticity and lower speed, high-angle flight conditions. Ad- 
vances have paralleled developments in steady Computational 
Fluid Dynamics (CFD) with a lag of approximately five years 
due to the additional requirement of time-accuracy. This pa- 
per presents a discussion of current aeroelasticity problem areas 
or challenges. The focus is primarily upon methods aimed at 
the study of nonlinear fluid dynamic flows, typically referred to 
as Computational Fluid Dynamics (CFD), although attention is 
also given to linear flow models. 

Figure 1 (Edwards [53]) illustrates significant features which 
must be addressed in the use of computational aeroelasticity 
for flutter boundary prediction. In this figure, a typical flut- 
ter boundary curve, characterized by the flutter speed gradually 
dropping to a minimum in the transonic speed range followed 
by a rapid upward rise, is shown. The ability to predict this min- 
imum, termed the transonic flutter dip, is of great importance in 
design, since the flutter boundary must be shown by a combina- 
tion of analysis and flight test to be outside the flight envelope 
by a specified margin. For military aircraft, the margin in terms 
of equivalent airspeed is at least 15 percent. Subsonic linear 
unsteady aerodynamic theories have been reasonably success- 
ful in predicting this flutter boundary for Mach numbers up to 
0.6-0.7 but linear theory is unable to account for the effects of 


aerodynamic shape and maneuvering condition upon unsteady 
airloads. At high Mach numbers linear analysis has been used 
with more or less success depending upon the severity of local 
transonic effects. The occurrence of flutter within the flight en- 
velope of an aircraft usually leads to structural failure and loss of 
the vehicle, highlighting the necessity for careful validation of 
computational methods intended for use in this area. In addition, 
aircraft service life can be significantly degraded by unforeseen 
dynamic loadings, such as buffet, and predictive capability for 
such off-design point loadings must be well-understood before 
being utilized in structural design. These key differences in the 
utilization of steady and unsteady computational methods must 
be clearly understood. 

This field received an initial impetus in the mid-1970’s from 
three sources: Tijdeman’s [152] experimental work on transonic 
unsteady pressure measurements, Magnus and Yoshihara’s [108] 
demonstration of key transonic flow features for an airfoil with 
an oscillating flap and the introduction of an economical tran- 
sonic finite-difference solution algorithm (LTRAN2) by Ballhaus 
and Gooijian [17]. Ballhaus [16] gives a survey of the field from 
this period. The AGARD Structures and Materials Panel Sub- 
committee on Aeroelasticity has selected experimental unsteady 
pressure data sets and defined two- and three-dimensional Stan- 
dard Aeroelastic Configurations [30, 31] to provide reference 
computational test cases for the development and validation of 
improved computational methods. 

Unsteady aerodynamics has been the theme of six recent 
conferences [8-12, 153] whose proceedings contain a wealth of 
information. Summary papers of the 1984 and 1985 AGARD 
conferences are given by Mykytow [114] and Mabey and Cham- 
bers [105]. The latter reference makes recommendations re- 
garding computational and experimental methods for unsteady 
flow phenomena and draws particular attention to the need to 
pay careful attention to the nature of shock motions. The peri- 
odic oscillations about circular arc airfoils are recommended as 
benchmark computational cases for all time-dependent transonic 
viscous flow theories. Zwaan [167] surveys aeroelastic problems 
in transonic flow while Deiwert [46] reviews the numerical sim- 
ulation of unsteady interactive flows. Reference [155] provides 
a collection of articles going into extensive detail for unsteady 
transonic aerodynamics. Mabey [106] gives a review of per- 
tinent experimental research on time-dependent aerodynamics. 
Finally, Dowell [48] provides an overview of nonlinear aeroe- 
lasticity phenomena including structural as well as aerodynamic 
nonlinearities. 
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Figure 1 Graphical Representation of Minimum 

Required Flutter Margin for Military Aircraft [7]. 


Figure 2 Characteristics of Attached and Separated 
Flow for Complete Aircraft (Edwards [53]). 


Research in these areas requires the comparison of exper- 
imental and computational results with the goal of achieving 
accurate predictive capability. Edwards [52, 54] provides sur- 
veys of these efforts for the transonic flutter problem while 
Mabey[107] discusses the physical phenomena associated with 
unsteady transonic flow. Bobbitt s [36] review of the issues in- 
volved in obtaining accurate results from experiment and from 
computation is particularly noted. Regarding higher angle, vor- 
tex dominated flows, a trend of increasing interest by the aerody- 
namics community in unsteady flows is also noted. This is due 
to the inherent unsteadiness of such flows and to the ability of 
emerging CFD methods to simulate their details. Newsome and 
Kandil [117] discuss physical modeling issues involved in the 
computational prediction of vortex dominated flows and survey 
numerical results. 

The remainder of this paper will review the current status of 
computational methods for unsteady aerodynamics and aeroelas- 
ticity. The key features of challenging aeroelastic applications 
are discussed in terms of the flowfield state: low-angle high 
speed flows and high-angle vortex-dominated flows. Next the 
computational methods and the basic fluid dynamic equations 
are introduced, followed by an assessment of the computer re- 
sources required for the unsteady aerodynamic computations. 
Then, the current state-of-the-art in CFD methods for transonic 
flows and vortex-dominated flows are each discussed, with em- 
phasis in the progress achieved during the past half decade. Fi- 
nally, an assessment of current capabilities and future research 
trends is offered. 

2. Features of Low-Angle, High Speed Flows 

The main features of steady transonic flow are described 
first in order to organize the discussion. With increasing Mach 
number and moderate angle-of-attack, the flow over the upper 
surface of an airfoil becomes critical between A/oo = 0.4-O.7 
with the first shocks forming at approximately 0.1 higher Mach 
number. Pearcy et ai. [120] have classified several types of 
flow separation which may occur. For conventional airfoils 
the typical pattern involves the growth of a local separation 
bubble, induced by boundary layer separation at the shock 
foot, spreading rapidly to the trailing edge as Mach numbers 
increases. This condition is often accompanied by unsteady 
phenomena such as buffet and aileron buzz (Tijdeman [152]). 


The steep aft pressure gradients of modem airfoils can lead 
to an alternate pattern in which separation progresses from the 
trailing edge towards the shock. Tijdeman [152] notes that 
the flow conditions in the region between the onset of trailing 
edge separation and fully separated flow are very sensitive to 
Reynolds number and the location of transition from laminar to 
turbulent flow. 

Figure 2 shows a diagram, from Edwards [53], of attached, 
mixed and separated flow regions for a complete aircraft at 
freestream Mach numbers between 0 and 2.0. In region 1, the 
flow is predominantly attached. To obtain optimum performance 
and to avoid the drag penalty associated with flow separation, 
design cruise conditions for aircraft typically are located here 
near the boundary of region II. 

As speed and/or angle of attack increase, a transition re- 
gion of mixed flow (region II of fig. 2) is encountered. For 
rigid structures, this region is typified by the onset of local- 
ized regions of flow separation which may exhibit significant 
aerodynamic unsteadiness. For realistic flexible structures, the 
aeroelastic response of the structure interacts with the airflow to 
induce much more complicated situations. For instance, struc- 
tural vibrations can cause the flow to alternately separate and 
reattach at flow conditions where a rigid structure would support 
attached flow. The associated high unsteady aerodynamic load- 
ing can interact with the structure to cause unusual aeroelastic 
phenomena which may restrict the vehicle flight envelope. 

With further speed and/or angle of attack increases which 
may be encountered under maneuvering conditions, fully sep- 
arated flow conditions emerge (region ID of fig. 2). Leading- 
edge vortex flows and stalled wing flows are of this nature. At 
still higher angles, vortex bursting in the vicinity of the aircraft 
can cause severe buffeting. Within such regions the flow is 
highly unsteady and accurate computations will require careful 
attention to turbulence modeling. To emphasize the complexity 
which the aeroelastic response adds, the flow within the three 
regions of figure 2 will be referred to as Type I, II, and III 
respectively. 

While the predictive methods for attached flows are reason- 
ably well developed, the picket fence in figure 2 emphasizes the 
difficulty in predicting aeroelastic phenomena in the mixed and 
separated flow regions. It also symbolizes novel features that 
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arc being encountered in transonic flutter testing. Modem high 
performance aircraft are capable of maneuvering at transonic 
speeds, leading to a much enlarged parameter space that must 
be considered in flutter analysis and testing. Wing/store loading, 
fuselage interference, angle-of-attack, Reynolds number, wing 
shape, and wing sweep all must be considered, and the tradi- 
tional flutter boundary parameterization of dynamic pressure at 
flutter versus Mach number may need to be augmented to ad- 
equately describe aeroelastic stability boundaries. For instance, 
flutter tests give some indication that these additional param- 
eters affect the detailed aeroelastic stability condition near the 
flutter boundary. Thus, the pickets of the fence in figure 2 rep- 
resent possible regions of low damping or instability that might 
be encountered 


Flutter 

dynamic 

pressure 



3. Features of High-Angle, Vortex — Dominated Flows 

Unsteady airloads due to flow separation are involved in a 
number of cases critical to the structural integrity of aircraft. As 
speed increases for moderate angles of attack, typical of maneu- 
vering flight near trimmed flight conditions, local transonic flow 
effects are encountered which lead to separated flow over the 
aft portions of lifting surfaces. Minimum flutter speed indices 
are often encountered in this transonic region, in conjunction 
with the onset of separated flow. The ability to predict these 
minimums is obviously tied to the ability to treat such “local” 
separated flows on wings. 

For slightly lower speeds where more aggressive maneu- 
vering is possible, unsteady airloads due to flow separation over 
“remote” components (e.g. forebody and main wing panel) leads 
to issues of tail buffeting and structural fatigue of aircraft com- 
ponents. For these cases, as speed, and/or angle of attack in- 
crease, smooth air flow over lifting surfaces breaks down in a 
variety of ways depending strongly upon the geometry. For 
lower sweep angles and blunt leading edges, flow separation 
may initiate near the trailing edge or near shocks and progress 
to completely separated and stalled conditions. For higher sweep 
angles and less blunt leading edges, leading edge flow separa- 
tion bubbles foreshadow the development of leading edge vortex 
flows. At higher angles, unsteady and burst vortex flow in the 
vicinity of the wing and downstream lifting surfaces leads to 
strong unsteady airloads and buffeting. Flow conditions near 
the boundaries of these regions for the various flow phenomena 
can be sensitive to a number of conditions and an understanding 
of these effects is called for in order to avoid adverse aeroelastic 
effects such as stall flutter, buzz, and structural buffeting. 

4. Computational Aeroelasticity Challenges 

With figure 2 providing a framework within which typical 
flowfields encountered in aeroelasticity may be viewed, a num- 
ber of current aeroelasticity problem areas are introduced in 
Table I and figures 3-7. Table I serves as a guide for discussing 
the current status in this area and the likely future trends. On 
the left are listed the key Challenges, most of which have been 
extensively commented on above. They are roughly graded in 
terms of increasing difficulty from top to bottom with the more 
difficult areas calling for more sophisticated flowfield modeling 
in order to achieve useful accuracies. Arrayed against these 
challenges are the Resource Issues impacting the economics of 
aeroelastic analysis, which are listed on the right. The choice of 
the appropriate level of CFD code to use, indeed the decision of 
whether to use a linear or nonlinear flow method, is dictated by 


Figure 3 Features of Transonic Flutter (Edwards [53]). 


the stage of design maturity, the required level of accuracy and 
computer resources available. The process of identifying criti- 
cal loading conditions requires running large numbers of cases, 
almost always utilizing lower level methods. Critical cases, so 
identified, are then candidates for further analysis with higher 
level methods. It is interesting to query if this process can be 
relied upon to capture the actual critical loading cases. 

Figure 3, from Edwards [53], indicates further features 
of high speed, low angle flutter. Dynamic pressure at flutter 
tends to decrease with increasing Mach number to a minimum 
“critical flutter point” value in the transonic speed range. At 
subsonic speeds where the flow can be assumed to be attached 
(Type I flow) at flutter, linear theory is reasonably accurate. 
As speed increases into the transonic region, the situation is 
complicated by the formation of shock waves and the onset 
of flow separation (Type II flow) and linear theory must be 
used with caution. The tow damping region in the figure 
indicates the potential for nonclassical aeroelastic response and 
instabilities which may be encountered. Figures 4-b, illustrate 
several types of novel aeroelastic responses which have been 
encountered with the onset of Type II flows and which offer 
challenges for computational methods. Figure 4 (Edwards [53]) 
shows a region of nonclassical aeroelastic response observed 
on a high aspect ratio, flexible, supercritical wing (Seidel et al. 
[141]) where high dynamic response at nearly constant Mach 
number was encountered at dynamic pressures well below those 
for which flutter was predicted with linear theory. The motion 
is of the limit-amplitude type and the response is believed to 
be associated with flow separation and reattachment driving the 
wing motion in the first bending mode. Figure 5 (Edwards [53]) 
illustrates wing/store limited amplitude oscillations experienced 
by modem, high performance aircraft under various loading 
and maneuvering conditions at transonic Mach numbers. Such 
oscillations can result in limitations on vehicle performance. 
The conditions for which this type of response occurs appear 
to also be near the onset of Type II mixed flow. The response 
typically increases for maneuvering flight conditions. Dynamic 
vortex- structure interactions causing wing oscillations have been 
observed, figure 6 (Dobbs et al. [47]), on a bomber type aircraft 
for high wing sweep conditions during wind-up turn maneuvers. 
The flow involves the interaction of the wing vortex system 
with the wing first bending mode and occurs over a wide Mach 
number range at moderate angles of attack. 
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For lower speed flight where higher angles are achieved, 
fully separated flows (Type III flows in fig. 2) are encountered 
which can range from diffuse vortical flow structures to concen- 
trated vortices designed to enhance stability and control. Inter- 
action of such forebody and wing vortex systems with aft ve- 
hicle components results in vortex-induced buffet loads. Figure 
7 (Edwards [53]) shows typical operating conditions at which 
such empennage buffet may be encountered. Buffet of horizon- 
tal tails can occur at intermediate angles of attack and is a result 
of the vortex system propagating downstream and encountering 
the horizontal tail surface. As angle of attack increases, the 
location of vortex bursting moves upstream in the wake. Loss 
of lift is associated with the burst location reaching the vicinity 
of the aircraft, and vertical surfaces located in such regions can 
experience severe dynamic loads and structural fatigue. 

These challenges, illustrated by figures 3-7, involve two 
types of unsteady flows. The first is the Type II flow (fig. 
2) wherein the onset of flow separation at high speeds leads 
to critical flutter conditions and/or novel aeroelastic responses. 
The second involves fully separated Type III vortex-dominated 
flows at high angles. The search for the appropriate levels of 
sophistication in fluid dynamic modeling to adequately model 
these flows is the subject of this paper. 


Table 1 COMPUTATIONAL AEROELASTICITY 
CHALLENGES AND RESOURCE ISSUES 


Challenges 

Resource Issues 

Stability & Control 

Modeling Tradeoffs 

divergence 

higher level CFD 

roll performance 

improved configuration 

wing rock 

detail 

Gust Response 

Design Maturity 
preliminary design 

Flutter Boundary Prediction 
l‘g 

maneuvering 
limit cycle oscillations 

final design 

critical loading conditions 
Computer Resources 

Control Effectiveness 

required level of accuracy 

buzz 

cost per solution 

hinge moment 

Buffet Response 
local: main wing panel 
remote: tail buffet 

number of solutions 


Dynamic 

pressure 

psf 



Figure 4 Region of High Dynamic Response Encountered 
During Test of a Flexible Supercritical 
Wing (Edwards [53]). 




Mach number 


Figure 5 Region of Potential Limited-Amplitude Oscillations 
Due to Wing-Store Loading Effects (Edwards [53]). 


Lift 

Coefficient 



Figure 6 Vortex-Induced Wing Oscillations of a Bomber-Type 
Aircraft (Dobbs, Miller, and Stevenson [47]). 


Figure 7 Regions of Vortex-Induced Buffet 
Loads (Edwards [53]). 
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5. Computational Methods 

A variety of fluid dynamic flow models is available to ad- 
dress unsteady aerodynamic computations. The choice of an 
appropriate method calls for assessment of the difficulty of the 
aerodynamic problem being addressed Type I flows (fig. 1) in- 
clude one of the most important aeroelastic analysis conditions, 
that of cruise at high dynamic pressure. Classical linear aeroe- 
lastic analysis has been primarily focused upon this condition. 
The transition from Type I to Type II conditions can occur due to 
aircraft maneuvering, with little decrease in dynamic pressure. 
Thus, aeroelastic response and stability of aircraft operating in 
Type II flows can be quite important although they have only 
begun to be brought within the range of computational methods. 

Fluid dynamic flow models available for unsteady aerody- 
namic computation include: the classical (linear) small distur- 
bance potential equation, the nonlinear potential equation (both 
Transonic Small Disturbance (TSD) and Full Potential equation 
(FPE)), the Euler equations (EE) and the Navier-Stokes equa- 
tions (both Full (FNS) and thin-layer (TLNS)). 

Issues which have been central to unsteady CFD have been 
the choice of implicit versus explicit algorithms, the stability 
of alternative solution algorithms and the treatment of compu- 
tational grids. Explicit schemes are simple to code and easily 
vectorizable but are limited in allowable time step by the sta- 
bility limit imposed by the signal propagation time over the 
smallest grid cell. Faced with the requirement of maintaining 
time-accuracy throughout the entire field for aeroelastic compu- 
tations, this easily leads to excessive computation times, espe- 
cially for viscous flow calculations where a very fine mesh near 
the surface is required to resolve the boundary layer. The al- 
ternative implicit solution algorithms thus are currently favored 
for present-day computer architectures because of their relative 
stability and time-step characteristics. While no attempt will be 
made to present complete details of the various levels of flow 
models, the following sections highlight the key equations and 
relevant boundary conditions. 

5.1 Navier-Stokes Equations 

Anderson, Tannehill and Pletcher [2] provide a description 
of the three-dimensional Full Navier-Stokes equations. For en- 
gineering applications, the Reynolds- averaged form of the FNS 
equations are normally used as a basis for practical, computa- 
tional procedures. In a cartesian coordinate system, the FNS 
equations can be written as follows: 

where the vector of independent, conserved variables is 

Q = (p,pu,pi\pw,e) T (2) 

and T 

F = ( pu , pu 2 + p, pur, puir, u(e + p)) (3) 

G = (pr, pur, pr 2 + p, pvxv, r(e + p)) (4) 

T 

II = (pw, puw, pvw, pw 2 + p, ir(e 4- p)) 


are the inviscid flux vectors. The viscous flux vectors, F V1 G v 
and H v are given in Ref. [2]. 

To facilitate solutions of the FNS equations on body-fitted 
computational grids, a curvilinear coordinate system is defined 
and eqs. (1) are converted using a generalized coordinate 
transformation of the type: 

{ = V = <T = c(x,y,r,/), T = t 

( 6 ) 

into the following expressions: 

l(«) + %(/■■ - f -) + r,( e - 1 *•) + M* - *•) = 0 

(7) 

where: 

Q = j(p,pu,pv,pw,e) T ( 8 ) 

F= ^j(pU,pUu + £xp, pUv + £ y p,pUw + U(e + p) - itp) T 

(9) 

1 T 

G = -{pV,pVu + T] z p,pVv + r}yp,pVw + rj 2 p,V(e + p) - rpp) 

J (10) 

h = lj{p\v,p\Vu + (; x p,p\Vv + (;yp,p\Vw + (;;p>\v(e + p)-Qp) T 

J (ii) 

Similar transformations are applied to the viscous flux vectors. 

Note that the transformation is time-dependent, allowing the grid 

to move to follow body motion and giving rise to grid motion 
terms such as £/, rp and ( t in eqs. (9) to (11). For viscous- 
flow aerodynamic computations, the solid-surface boundaries 
are modeled using the “no-slip” condition, together with adia- 
batic wall and zero normal pressure gradient conditions. 

A modified form of the FNS equations, termed the thin- 
layer Navier-Stokes equations, has been found useful for ap- 
plications where viscous effects in certain spatial directions are 
small enough to be neglected. For many aerodynamic flows of 
interest, the viscous terms normal to the body are of most im- 
portance, and the other viscous fluxes can be dropped. If the 
^-direction in eq. (7) is taken as the body normal direction, a 
TLNS form of eq. (7) is given by the expression: 

l (^) + ^( / ’) + ^( 6_6 ”) + ^( // ) =0 (12) 

For turbulent-flow calculations, turbulence modeling such as 
the algebraic eddy viscosity model of Baldwin and Lomax [15] 
is used. Rumsey and Anderson [134] are typical of applications 
using this thin-layer approximation to compute viscous-flow 
solutions for airfoils. Also, Thomas et al. [150] describe a 
three-dimensional implementation of the above equations in the 
CFL3D code. 

5.2 Euler Equations 

For aerodynamic flows in which viscous effects are expected 
to be negligible, the inviscid Euler equations can be derived 
from eq. (7) by dropping all three of the viscous flux vectors 
from the formulation. Then, F v = G v = H v - 0, and eq. (7) 
reduces to: 

l («) + |( # ) + i |( i ) + ;|(*)= 0 <13 > 


(5) 
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The boundary conditions appropriate to the Euler equations 
are the “slip” or “flow tangency” conditions. For these applica- 
tions, only the velocity component normal to the body surface 
is set to zero. The flow streamlines are assumed to run parallel 
to the surface tangent at each point on the surface. Note that 
for situations where there is rotational flow, such as the regions 
behind strong or curved shock eaves, the Euler equations can 
propagate the vorticity downstream in a correct manner if an 
adequate number of grid points are used. 

5.3 Full Potential Equation 

The FPE is derived from the Euler equations by assuming 
that the flow is inviscid, isentropic and irrotational. A velocity 
potential can then be defined whose derivatives in the spatial 
directions recover the flow velocity components in the appro- 
priate directions. The governing equations for this formulation 
are the continuity equation: 

Pi + (p*x)x + (/»*,), + (/»*,), = 0 (14) 

and the isentropic energy equation: 

(15) 


l + V^'&O 


24-, *;) 


As with the Euler equations, a generalized coordinate system is 
often used to solve the FPE and a flow-tangency boundary con- 
dition is enforced at the body surface. Sankar and Malone [136] 
describe such a FPE formulation in generalized coordinates. 


wake. The pressure coefficient may be computed using either 
linear or nonlinear forms of the Bernoulli equation. Batina et 
al. [21] describe this algorithm as implemented in the CAP-TSD 
code with a number of example calculations. 

5.5 Viscous -Inviscid Interaction 


Neither the potential equations nor the Euler equations de- 
scribed above incorporate viscous effects which can be impor- 
tant for high speeds and for lower speed at higher angles. It is 
possible to account for unsteady viscous effects by coupling a 
viscous boundary-layer model with an otherwise inviscid anal- 
ysis. As commonly implemented, the inviscid outer flow solu- 
tion provides the surface pressure distribution needed to solve 
the boundary layer equations. This yields the boundary- layer 
displacement thickness distribution which is used to modify the 
airfoil surface tangency boundary condition for the next iteration 
of the outer inviscid flow solution. 

Guruswamy and Goorjian [65], Howlett and Bland [81], 
and Rizzetta [129] describe this method implemented in two- 
dimensional unsteady TSD codes. The effect of a viscous 
boundary layer for attached turbulent flow is modeled in a quasi- 
steady manner by means of the lag-entrainment equations of 
Green et al. [64]. In this integral method the displacement 
thickness 6* is computed as a function of the boundary -layer 
momentum thickness 6 and the shape factor H as 

6'=0H ( 21 ) 


5.4 Transonic Small Disturbance Potential Equation 


The Transonic Small Disturbance Potential equation is de- 
rived from the inviscid Euler equations assuming that the flow 
is isentropic, irrotational and a small perturbation of a steady 
uniform flow, U , in the x direction. The TSD velocity potential 
function, <5, describes the perturbed velocity components u, v, 

d$ d& do 

u = — , v = — , u; = — 

ox ay dz 


(16) 


where the total velocity in the x direction is U+u. Batina [20] 
and Batina et al. [21] give the modified TSD potential equation 
in conservation form as 


Ofo dh dh d_h_ 
di*dx*dy dz 

where 


(17) 


/o - —A&t - 0$ x ; fo = 4> y + H<t> z Q y 

fi = Eo z + F6i + G6] ; /, = <j> z 

(18) 

The coefficients A-H are given by Batina [20]. The TSD 
equation (17) is distinguished from the higher equation level 
flow models in that, within the small disturbance assumption, the 
computational grid is not required to move with the body since 
boundary conditions are imposed at the mean plane, usually z 
= 0*. The wing flow tangency boundary condition is 


= ft + ft 09) 

where / i (x,y,t) = 0 describes the upper and lower body sur- 
faces. The trailing wake boundary conditions are 


(tfx + &] = 0 ; [<t> 2 ] = 0 (20) 


where [-] indicates that jump in the indicated quantity across the 


Given the velocity at the edge of the boundary layer (from the 
outer inviscid flow solution), the boundary layer equations may 
be integrated in a “direct” fashion to obtain 6*. Coupling be- 
tween the boundary- layer and the outer inviscid flow is through 
the boundary conditions on the airfoil and wake, eqs. (19) and 
(20), which are modified to 

= ft + ft + (6'/c)x ; M = P7c)J (22) 

For cases of incipiently separating and separated flows the 
boundary layer equations become singular, requiring a refor- 
mulation of the equations in an “inverse” mode in which the 
edge velocity gradient is computed for a given displacement 
thickness (East et al. [49]). Consistency with the outer inviscid 
flow may be obtained via a “semi-inverse” relaxation coupling 
method (in which 8* is updated based upon the error between 
inner and outer edge velocities) described by Carter [38] and 
also used by Fenno et al. [58]. For cases with large amounts 
of separated flow, particularly for unsteady flow, the semi- 
inverse method itself encounters stability problems (Edwards 
and Carter [50]). These cases have been more tractable via 
the “quasi-simultaneous” coupling method introduced by Veld- 
man [158] and by Houwink and Veldman [79] and the “semi- 
implicit” coupling method of LeBalleur and Girodroux-Lavigne 
[98]. Both of these methods perform the viscous-inviscid cou- 
pling by developing, at each grid point, locally linear relations 
between the inner and outer flow variables. This enables si- 
multaneous solution for the coupling variables which is usually 
accompanied by relaxation and iterations for convergence. The 
quasi-simultaneous method has been implemented using the low 
frequency LTRAN2-NLR TSD code and quasi-steady integral 
boundary layer equations. The semi-implicit method described 
in Ref. [98] achieves full time-consistency by coupling a time- 
accurate TSD code with a time-accurate integral boundary layer 
method. 
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5.6 Time-Linearized Transonic Small Disturbance Equation 

A second formulation of the Transonic Small Disturbance 
Potential equation is the time-linearized equation, which is de- 
rived by assuming that the unsteady flowfield can be treated as a 
small perturbation about the steady flow field solution. This the- 
ory assumes that shockwaves are neither created nor destroyed 
during the unsteady motion. The steady flow potential is ob- 
tained from solutions to the steady-flow version of the TSD 
equation: 

[1 - Ml - (7 + 1).VM + t° yy + = 0 (23) 

The unsteady potential, is then computed from the unsteady 
TSD equation by solving the following equation: 

^i,+{ [1 - Mi - (7 - 1 = 0 

( 24 ) 

Note that the steady potential <^°, is required in the above 
equation. In practice, <5° can be obtained from other theoretical 
formulations or derived from experimental data. Hounjet [74] 
is representative of this approach to unsteady transonic flow 
modeling. 

5.7 Computational Grid Effects 

Grid generation for unsteady problems in which the body 
boundary moves, such as for an oscillating control surface 
or an aeroelastic deformation, raises new issues over those 
involved in steady flows. At the TSD equation level, wherein 
the boundary conditions are applied on a nonmoving, mean 
surface plane, stationary Cartesian grids have been used. For 
unsteady problems, care must be taken to ensure the fidelity 
of signals propagated through the stretched grids which are 
used. Seidel et al. [140] and Bland [34] give results for the 
TSD equation. A key effect of grid stretching is to modify 
the “dynamic impedance” of the mesh at internal grid points, 
leading to “internal reflections” of waves which may return to 
the vicinity of the modeled aircraft components and contaminate 
unsteady solutions. It is shown that this issue is typically more 
severe for exponentially stretched meshes, frequently used for 
steady calculations, than for meshes with less severe stretching 
in the near field. This effect is alleviated in three-dimensional 
flows. Bland [35] provides guidelines for generating grids and 
selecting time-step size for accurate unsteady computational 
characteristics. Finally, this effect has not been observed or 
documented for calculations utilizing the higher level Euler or 
Navier-Stokes equations. 

For flow modeling equations higher than the TSD equation, 
the body-conforming grids used must be realigned with the mov- 
ing body at each time step to maintain accuracy. Schemes for 
accomplishing this have been studied as well as the necessity 
of moving the grid at all. When body motions are small with 
perturbations mainly normal to the surface, imposing “transpira- 
tion” boundary conditions on the mean surface location may be 
an acceptable approximation (Sankar et al. [137]). Steger [147] 
formulated the TLNS equations including terms accounting for 
grid motion. Steger and Bailey [148] used simple shearing of the 
grid coordinate normal to the surface to allow the grid to follow 
aileron motions. Chyu and his coworkers [42, 43] used an inter- 
polation scheme for defining grids at intermediate steps between 
the extremes of motion for oscillating airfoils and for fixed outer 
computational boundaries. Anderson et al. [3] present EE re- 
sults for dynamically moving airfoils and wings in which the 


entire mesh is rotated to follow rigid airfoil and wing motions. 
For aeroelastic motions of flexible structures more general meth- 
ods for dynamically moving the mesh are required. Guruswamy 
[67], Ide and Shankar [84], and Nakamichi [115] describe meth- 
ods wherein the curvilinear coordinate normal to the surface is 
sheared based upon the instantaneous surface normal displace- 
ment as computed by the time-marching aeroelastic equations. 

The above methods have all been implemented using struc- 
tured grid meshes in which computed variables for neighbor- 
ing grid points are stored in adjacent computer memory cells. 
Unstructured grids, which can be implemented with triangular 
grid cells in two-dimensions and as tetrahedral cells in three- 
dimensions, offer more flexibility in modeling complex geome- 
tries. Batina has developed a method for moving such body 
conforming meshes to maintain alignment during aeroelastic 
motions for airfoils [24] and complete configurations [23]. A 
network of springs is associated with the mesh in which the 
edge of each mesh cell is represented by a spring whose stiff- 
ness is related to the length of the edge. At each time step, the 
new location of the body boundary obtained from the aeroelastic 
equations of motion is used to solve for the new static equilib- 
rium location of the nodes of the spring network. Robinson et al. 
[130] modified this grid motion technique for structured grids 
and give results of wing flutter calculations using an Euler code. 
Rausch et al, [124] further refined the method, treating spatial 
grid cell adaption (mesh enriching and coarsening procedures) 
on unstructured meshes. 

6. Computational Aeroelasticity Resource Assessment 

In this section, estimates of computer resources necessary 
to produce accurate, converged results are given. This is done 
prior to the discussion of the detailed applications in the fol- 
lowing sections in order to provide a framework within which 
to assess what has been accomplished and where further work 
is necessary. 

The measure which will be used for computational aeroe- 
lasticity resource requirements is the Computer Processing Unit 
(CPU) run time for a single steady or dynamic time-marching 
calculation. Typically, a steady calculation is required to es- 
tablish an initial flowfield for subsequent dynamic calculations. 
Each dynamic calculation produces a time history of aeroelas- 
tic response from which stability or response measures, such as 
modal frequency and damping, may be derived. The CPU run 
time for a CFD calculation can be estimated from the relation 

Tcpu = N C yN s tf C yNgpT (25) 

relating computer CPU run time, T cpu , to the number of com- 
putational steps, N a t\ the number of cycles of oscillation for a 
given frequency, N cy \ the number of steps per cycle (required 

for accuracy or stability), N at / Cy \ the number of grid points, 
N gp \ and the algorithm speed, r. 

The algorithm speed, r, is a common measure of the speed 
of an algorithm given in microseconds per grid point per time 
step (/isec/gp/st). Values used for this parameter herein assume 
machine speeds typical of a Cray 2 class supercomputer: 250 
million floating point operations per second. Lower values of 
r are associated with less complex algorithms, such as explicit 
methods, while more complex algorithms yield larger values. 
However, the higher level algorithms (e. g. implicit, upwind- 
biased, etc.) allow larger time steps and are generally favored 
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Table II COMPUTATIONAL RESOURCE ESTIMATES FOR HIGH SPEED, LOW ANGLE CONDITIONS 
Single half-span wing; up to “locally separating” flow 

T C pu — A eyA st/cyNgpT 




TSD 

TSD-VISC 

Euler 

TLNS 

Algorithm speed. 

T 

3.0 

6.2 

60 

65 

/zsec/gp/st 






Number of steps 






Steady state 

N st 

1000 

1000 

1000 

5000 

Dynamic 

Ncy 

3-7 

3-7 

3-7 

3-7 


N st /cy 

100-300 

100-300 


1000 

Number of grid pts 

N gP 

200 K 

400 K 

240 K 

650 K 

Run time, hrs 
Steady 

Tqju 

0.2 

0.7 

4.0 

59* 

Dynamic 


0.05 

0.2-1.4 

■ 

3. 6-8.4 

35-82* 


* Assumes uniform global time step; no zonal decomposition 


for calculations requiring time -accuracy. For a TSD code, 
values are in the range of 3-10 //sec/gp/st while Navier-Stokes 
codes call for values ranging from 10-100 //sec/gp/st. A value 
of 65 will be used below for the Navier-Stokes code resource 
estimates. This is representative of the speed for a three- 
dimensional implicit, upwind-biased code. 

Table II presents computational resource estimates for a sin- 
gle steady or unsteady (aeroelastic) response calculation for a 
high speed, low-angle problem. This implies that attached or 
only mildly separated flows are to be treated. The estimates are 
for a single half-span wing panel which is typically utilized when 
symmetry may be assumed. The four classes of CFD codes 
listed are: inviscid TSD, TSD with interacted boundary layer 
model (TSD-VISC), Euler equations, and Thin-Layer Navier- 
Stokes equations. The number of steps for steady and dynamic 
response calculations are typical of those reported in the litera- 
ture. The number of grid points have been selected based upon 
published grid convergence studies: see Anderson and Batina 
[4] for TSD and EE results and Vatsa and Wedan [157] for Eu- 
ler and Navier-Stokes results. The second and fourth columns 
giving estimates including viscous modeling are the most in- 
teresting as this level of flow modeling is required in order to 
achieve desired levels of improvement over well developed lin- 
ear unsteady aerodynamic methods. The estimates range from 
approximately one hour per case for the TSD-VISC capability 
to over 50 hours for codes based on the TLNS equations. Treat- 
ment of full-span configurations would double these estimates 
and increased geometric complexity, such as additional lifting 
panels and bodies which would be called for in complete con- 
figuration modeling, would also increase these estimates. Each 

additional component modeled (tail, fuselage, etc.) would cost 
roughly the amounts given in the Table for the single wing. This 
last cost factor is due to the nature of grid point densities used, 
wherein the large majority of grid points are clustered near the 
body surfaces with the grid density coarsening quickly away 
from surfaces. It is noted that the TLNS estimate is likely to 
be conservative in that a uniform global time step (limited by 
small grid cells in boundary layers) is assumed. A number of 
straightforward modeling changes could lead to more econom- 
ical methods with no loss of accuracy. 


Higher-angle vortex-dominated flows call for significantly 
more complicated flow modeling: flow separation over the fore- 
body/wing must be accurately predicted; vortex formation and 
the convection of vorticity over significant distances without loss 
of accuracy are involved. Perhaps the most difficult computa- 
tional aeroelasticity challenge is vertical tail buffeting wherein 
vortical flow from the region of the forebody/wing is convected 
downstream, possibly encountering vortex bursting, and impacts 
on the flexible tail structure. CFD computations of these features 
are beginning to emerge and are reviewed below. Edwards [55, 
56] assesses the computer resource requirements for such cal- 
culations along the lines given above. The flowfield is assumed 
decomposed into three regions: an inner viscous region adjacent 
to bodies and lifting surfaces, a nearfield region encompassing 
the vehicle and the vortical flow region above it and an outer 
inviscid, irrotational region. Assuming these three regions to be 
modeled using TLNS, Euler and potential equations respectively 
leads to the following rough computer resource estimate for a 
tail buffeting calculation: 


T — N N , L maz - 

j-cpu — cy 1 ' alfcy r 1 
Jmin 


196,000 N bod + 


u* 


+ 50,000 
( 26 ) 


The frequency bandwidth of interest is specified by / min and 
/ mar, A <bod is the number of lifting surfaces and bodies wetted 
by viscous flow, Vjj is the volume of the nearfield vortical flow 
region and N gp f wl is the number of grid points per spatial wave- 
length assumed necessary for accurate vortical flow calculations. 
The sample estimate given by Edwards indicates run times on 
the order of 1000 hours for such a calculation. Use of an Euler 
code for a similar calculation leads to reduced estimates on the 
order of 100 hour run times. However, issues related to the 
use of inviscid Euler codes for this type of calculation indicate 
caution in their use [55]. 


7. Low-Angle, High-Speed Flow Applications 

This section will discuss applications of CFD methods for 
flows which are generally at lower angles and high speeds. 
First, the available experimental databases needed for valida- 
tion of CFD codes are described. Then, recent applications of 
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CFD methods for unsteady transonic aerodynamic analysis are 
presented, followed by a survey of aeroelastic applications for 
flutter prediction. Finally, the status of complete aircraft con- 
figuration modeling is illustrated through cited examples found 
in the current literature. 

7.1 Experimental Data Bases for CFD Code Validation 

The AGARD Structures and Materials Panel Subcommittee 
on Aeroelasticity has selected experimental unsteady pressure 
data sets as Standard Configurations to provide reference com- 
putational test cases for CFD code validation. Data from es- 
sentially rigid models oscillating in selected degrees of freedom 
is available for six airfoils (Bland [30]) and five wings (Bland 
[31]). Reference [11] contains many comparisons with these 
data sets. Edwards and Thomas [52] survey calculations for 
these and other cases, providing tables of example applications. 
Detailed comparisons of calculations with these experimental 
data sets can be found in Angelini et al. [5], Houwink [77], 
Bland and Seidel [32], and Howlett and Bland [81]. 

The major intended use of unsteady aerodynamic calcula- 
tions is for prediction of aeroelastic response of aircraft and, 
more specifically, flutter. There have been numerous published 
calculations of two degree-of-freedom airfoil studies which are 
devoid of experimental data comparisons since realistic 2-D 
flutter models are very difficult to fabricate. On the other hand, 
there are only a small number of published 3-D flutter cal- 
culations which are compared with experimental results. An 
important reason for this is the detail and effort required to per- 
form a valid flutter analysis of a flutter model. Vibration mode 
shapes and masses must be accurately calculated or measured 
and surface coordinates measured. 

In many instances, experimental flutter data obtained for 

actual or proposed flight vehicles are considered proprietary 

by individual private organizations, and consequently, are not 
in the public domain. However, Yates [166] describes an 
AGARD standard aeroelastic flutter model consisting of a 45 
degree sweptback wing for which extensive flutter test results are 
available. Also, NASA Langley Research Center has recently 
begun a multi-year experimental program to generate extensive 
flutter and unsteady aerodynamic data suitable for aeroelastic 
CFD code validation applications. This research activity, known 
as the Benchmark Models program, is outlined in Bennett et 
al. [27]. Rivera et al. [126] describe the first flutter and 
unsteady pressure measurement model test completed as part 
of this ambitious test program. 

7.2 Unsteady Aerodynamic Applications 


7.2.1 Linear Aerodynamic Methods. Although the major fo- 
cus of this survey is on computational methods for nonlinear 
unsteady flows, it is prudent to be aware of the capabilities 
of linear methods for two reasons. First, economy of effort 
demands that proof be offered that supposedly more accurate 
methods do indeed make a difference. Secondly, new methods 
should always be checked against older well calibrated methods 
for conditions at which both should produce similar results, thus 
protecting against inappropriate use of “higher-order” methods. 

For subsonic lifting surfaces, the standard linear methods 
are the doublet-lattice method (e.g. Rodden et al. [131] and 
the kemal function method (e.g. Rowe and Cunningham [133]). 


Recent developments of unsteady supersonic aerodynamic the- 
ories have been the potential gradient method of Hounjet [76] 
and the harmonic gradient method of Chen and Liu [40, 41]. 
These references contain numerous examples of aeroelastic ap- 
plications. 

These linear aerodynamic methods have been highly refined 
and calibrated with aeroelastic model test data. Coupled with 
powerful finite element structural dynamic analysis programs, 
they are used at all stages of design. Two notable analysis codes, 
NASTRAN and Elfini, are described by Johnson and Reymond 
[89] and Petiau and Brun [121]. 

Liu et al. [104] describe extensions to an unsteady aero- 
dynamics panel method to permit more accurate modeling of 
wing/fuselage configurations for subsonic flow field conditions. 

Improvements to the technique used to implement unsteady 
boundary conditions on the body surfaces, the addition of an 
embedded singularity to simulate the aft flow separation in the 
body wake region, and the use of constant-pressure wing panels 
(as opposed to doublet lattice panels) are discussed in the paper. 
Results are presented for wing-tail, body-alone, and wing/body 
configurations. 

7.2.2 Nonlinear Aerodynamic Methods. Much effort has been 
directed towards predicting unsteady transonic aerodynamic 
flows using the transonic small disturbance equation. The 
XTRAN3S TSD code has been extended to allow treatment 
of wing-tail (Batina [18]) and wing-fuselage (Batina [19] and 
Guruswamy and Goorjian [70]) configurations. Batina [20] has 
adapted an approximate factorization (AF) solution algorithm 
based upon the work of Shankar et al. [142] to the TSD poten- 
tial equation. The AF algorithm is implemented in a computer 
code termed CAP-TSD (Computational Aeroelasticity Program 
- Transonic Small Disturbance) developed at NASA Langley 
Research Center (Batina et al. [21]). The code permits the 
aeroelastic analysis of complete aircraft through the modeling 
of multiple lifting surfaces and bodies. Results are presented for 
five configurations illustrating this capability. Steady and un- 
steady pressures for the F-16C aircraft modeled by four lifting 
surfaces and two bodies are presented and agreement with exper- 
imental steady pressures is considered good. The grid used for 
these calculations contained 324,000 points. The calculations 
required 0.88 CPU seconds per time step or 2.7 microseconds 
per grid point per time step on the CDC VPS-32 computer. 
Thirteen million words of memory was required. Pitt et al. 
|123] give flutter analyses obtained with the XTRAN3S code 
and the CAP-TSD code. Both codes were used for the wing- 
alone analysis of the F-15 and F/A-18 aircraft. The latter code 
was used to study canard/wing/tail and wing/launcher/tip missile 
configurations. There is a general lack of unsteady experimen- 
tal data for complex configurations with which to validate such 
computations. 

Rodman, Nixon and Huttscll [132] describe modifications to 
the XTRAN3S code to permit the use of experimental steady- 
flow pressure data in the unsteady flow solution of the TSD 
equation. The TSD equation is split into a steady and an 
unsteady component Experimental surface pressure data are 
used to define a “strained-coordinate” system, which is then 
used to solve the perturbation equation for the unsteady flowfield 
potential. The method is limited to those situations in which 
shockwaves are preserved throughout the unsteady motion. Test 
cases are presented for both 2D and 3D configurations (NACA 
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64A006, NLR 2302, airfoils and the F-5 wing). 

Angelini, et al. [5] report on a large number of unsteady 
transonic-flow computations performed using five different CFD 
methods. The computational procedures include a 1) 2D TSD 
code, 2) a 2D coupled TSD code/boundary-layer method for un- 
separated turbulent flows, 3) a 2D coupled TSD/boundary-layer 
code for unseparated and separated laminar and turbulent flows, 
4) a 2D Euler code for inviscid flows, and 5) a 3D TSD code 
for inviscid flows. Computed results are presented for several 
AGARD standard airfoil configurations and a rectangular wing 
model. Some comparisons are made with available experimen- 
tal data. 

Sotomayer, Sankar, and Malone [146] provide a comparison 
of computed results for the F-5 wing using three different numer- 
ical procedures. Results are presented for a transonic small dis- 
turbance code (XTRAN3S), a full potential code (USIPWING), 
and an Euler code (GTEUL-3D). Steady and unsteady-flow re- 
sults are compared with experimental data for several Mach 
numbers (Afoo = .8, .9, .95, .99 and 1.35). Sotomayer and Bor- 
land [145] describe an application of the XTRAN3S code to 
the prediction of unsteady wing airloads arising from oscillating 
control surfaces. The boundary conditions used to simulate un- 
steady control surface motions are outlined. Numerical results 
for the F-5 wing with a deflected and oscillating trailing-edge 
flap are presented for a range of transonic Mach numbers. Com- 
parisons are also made with existing experimental data for this 
wing. 

A number of investigators have also used the FPE as a 
basis for developing unsteady transonic flow solution proce- 
dures. Isogai [86] describes a procedure to solve the quasi-linear 
form using a finite difference procedure on a stretched Cartesian 
mesh. A quasi-conservative form of Jameson’s rotated differ- 
encing scheme is used to capture shock waves. Since the grid 
used in this method is not boundary conforming, interpolation 
formulas are used to project values of the surface boundary 
conditions to appropriate computational grid points near the ac- 
tual body surface. A semi-implicit solution procedure is used 
to advance the resulting equations in time. Computed results 
and experimental data are presented for the ONERA M6 wing 
(steady flow) and the NORA wing (steady and unsteady flow). 

Sankar and Malone [136] describe two different procedures 
(SUNTANS and USIPWING) used to solve 2D and 3D prob- 
lems, respectively. These methods differ in the manner in which 
surface motion is prescribed. In the SUNTANS code [109] the 
computational grid deforms with time to permit an exact imple- 
mentation of numerical boundary conditions at the instantaneous 
airfoil surface. In the USIPWING code [135], a surface tran- 
spiration technique is used to simulate small unsteady motions 
of the wing surface. Reference [136] presents 2D and 3D com- 
puted results compared to experimental data for a number of 
airfoil and wing geometries. 

Additional 2D FPE methods are presented by Schippers 
[138] and Whitlow [163]. Schippers describes the mathematical 
formulation and numerical implementation found in the NLR 
TULIPS code. The TULIPS code solves the 2D full-potential 
equation in strong-conservation form. The FPE is solved on 
a body-fitted grid, using a second order time-integration pro- 
cedure, together with flux-splitting techniques. The resulting 
finite-difference equations are integrated in time using the im- 
plicit method of fractional-steps. The algorithm is demonstrated 


for unsteady transonic flows about the NACA 64 AO 10 airfoil 
and an airfoil of the ONERA M6 wing. Whitlow describes 
modifications to a 2D Full Potential equation solution procedure 
which simulate nonisentropic flow effects arising from strong 
transonic shockwaves. A nonisentropic formulation is used to 
modify the density upwind biasing incorporated in the original 
FPE code. Calculated results are presented for several airfoil 
shapes. 

The important effects of control surface effectiveness and 
hinge moments are studied in Bharadvaj [29] and Ominsky and 
Ide [1 19]. Both studies involve full potential equation codes and 
the former includes the effect of a 2-D quasi-steady boundary 
layer model to calculate pressures due to oscillatory control 
surfaces on a high aspect-ratio supercritical wing. The latter 
reference studies multiple surface control of a free-to-roll model 
for simultaneous flutter suppression and roll control. 

Damodaran [45], Whitfield et al. [162J, Chaderjian and 
Guruswamy [39], and Nakamichi [115] provide four exam- 
ples of the use of structured-grid Euler/Navier-Stokes solvers 
for unsteady aerodynamic applications. Damodaran describes 
an application of Jameson’s explicit Runge-Kutta time-stepping 
scheme to the solution of the 2D Euler equations for unsteady 
transonic flows about oscillating airfoil configurations. Nonre- 
flecting far-field boundary conditions are used to help reduce 
non-physical wave reflections at the outer boundary of the com- 
putational grid. Numerical results are presented for the NACA 
64 AO 10 and the NACA 0012 airfoils. Additional calculations 
are presented for an airfoil experiencing axial oscillations (simu- 
lating a helicopter-rotor environment) and for the transient aero- 
dynamic response of an airfoil due to angle of attack change. 
For the latter case, comparisons are made with classical lin- 
ear aerodynamic theory. In Ref. [162], the authors describe 
a solution procedure for the 3D Euler equations in transonic 
flow. Dynamic multi-blocked grids are used to model complex 
aircraft configurations. The Euler equations are cast in strong 
conservative form and solved using an implicit, approximately 
factored scheme. Computational results are presented for a 
wing with store and a wing/counter-rotating prop-fan configura- 
tion. Limited comparisons with experimental data are provided 
for the wing/store configuration. Reference [39] describes a 
computational procedure for the thin-layer Navier-Stokes equa- 
tions. The method can be used to model complex configu- 
rations via the multi-block technique for grid-generation. An 
implicit approximate-factorization diagonal algorithm is used to 
integrate the fluid dynamics equations. The method is con- 
sidered first-order accurate in time. Computational results are 
presented for a low-aspect-ratio rectangular wing. Both steady- 
flow and unsteady-flow calculations are compared, and one case 
is presented for which the wing grid is split into an inner and 
an outer region in order to demonstrate the time accuracy of 
the method. Finally, Nakamichi presents comparisons of TLNS 
code results for the NORA wing, one of the AGARD Standard 
Configurations. A moving grid capability is employed, allow- 
ing the computational grid to follow the harmonically pitching 
wing motion. 

Batina [23] has developed an Euler code based upon an 
unstructured grid with very general configuration modeling ca- 
pability. Results are given for a supersonic fighter configura- 
tion with canard/fuselage/cranked delta-wing/tail/flow-through 
nacelle components modeled. Calculations are given for the 
complete vehicle oscillating in an aeroelastic mode and utiliz- 
ing the moving mesh technique described above. Such unstruc- 
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tured grid methods involve more complex programming than 
structured grid methods. Efficiency can be regained by algo- 
rithm improvements such as the temporal-adaptive algorithm 
described by Kleb et al. [93] which may relax the stringent 
requirement of using a global time-step dictated by the size of 
the smallest grid cells. See also Angelini and Soize [6] for 
a two-dimensional small perturbation treatment of body-fitted 
unstructured grids for unsteady calculations. 

Mortchelewitcz and Sens [111] also describe a solution pro- 
cedure for the 3D Euler equations using unstructured grids. An 
explicit-implicit two phase integration scheme is used to solve 
the discretized equations on a mesh of tetrahedral cells. Surface 
transpiration is used to simulate motion of the solid surfaces for 
unsteady computations. Sample computations are presented for 
a wing-body configuration and compared to experimental data. 
The experimentally measured wing deformation is included in 
the computational mesh to account for static aeroelastic effects. 

7.2.3 Hybrid Aerodynamic Methods. A number of hybrid com- 
putational schemes have been developed which combine features 
of different aerodynamic theories. Liu and Kuo [102, 103] de- 
scribe a Transonic Equivalent Strip (TES) method for computing 
unsteady transonic flow about three-dimensional wing configu- 
rations. The method uses 2D steady-flow pressure distributions 
to determine an “equivalent” airfoil shape obtained by an auto- 
mated airfoil design procedure. This insures that the mean-flow 
transonic effects are matched closely. Then the equivalent shape 
is used to compute 2D unsteady transonic-flow using a nonlin- 
ear method, such as the LTRAN2 TSD code. Finally, the 2D 
unsteady transonic airloads are applied at each span station of 
the wing, using phase-lag modifications obtained from a 3D 
linear wave theory. Computed unsteady pressure distributions 
are given for several wing configurations, including wings with 
oscillating control surfaces. 

Hounjet [74] describes a hybrid field-panel/finite-difference 
procedure (FTRAN3) which combines features of both linear 
and nonlinear methods. The field panel procedure is used to 
compute a proper radiating boundary condition on the outer 
boundary of the grid used for the finite-difference method. This 
permits the grid extent to be reduced in size, resulting in a 
decrease in computational effort overall. Hounjet implements 
a multigrid technique to further improve convergence of the 
finite-difference method. 

Voss [160] describes other hybrid schemes. An unsteady 
field panel method (PTRAN3) for solving the time linearized 3D 
Transonic Small Disturbance equation is described. The wing 
and wake surfaces are modeled with unsteady dipole distribu- 
tions, while the region near the wing surface, where compress- 
ibility effects are greatest, is modeled using volume sources. 
Voss also describes a combined field-panel/finite-difference pro- 
cedure. 

Zimmerman and Vogel [169] describe a time-linearized 
method for solving the transonic small perturbation equation 
(LIN TSP). Comparisons are made for the unsteady transonic 
airloads computed using several different methods, including 
LIN TSP, TSP and doublet lattice. Some calculations are made 
using both a steady and a quasi-steady boundary layer. The pa- 
per also investigates the behavior of higher order harmonics of 
the computed airloads, and concludes that for small amplitude 
motions the higher-order harmonics can be neglected. 


Steiginga and Houwink [149] describe an engineering type 
method (Q3D) which couples 2D linear theory, 2D nonlinear 
theory and 3D linear theory, in order to predict unsteady tran- 
sonic airloads on transport-type aircraft wings. In this method, 
2D scaling matrices are first obtained by forming a ratio of un- 
steady 2D nonlinear to 2D linear AICs, where each nonlinear 
AIC is computed for a representative airfoil section of the 3D 
wing. The resulting scaling matrices are used to modify 3D 
AICs obtained from a doublet lattice analysis. 

Finally, Reisenthal and Nixon [125] have supplemented 
the TSD equation with a transport equation for streamwise 
vorticity and a vector potential equation to predict vortex effects 
over missile configurations. Results for a complete missile 
configuration at subsonic, transonic, and supersonic speeds are 
encouraging. 

7. 3 Aeroelastic Applications 

A restricted selection of the many publications of 2-D- 
potential flow flutter calculations includes Isogai [85], Edwards 
et al. [51], Bland and Edwards [33], Berry et al. [28], and 
Zimmerman [169]. References [33] and [51] document the 
somewhat surprising “locally linear” nature of transonic poten- 
tial flows. That is, about the steady mean flow condition which 
is a nonlinear function of Mach number and angle-of-attack, un- 
steady perturbation airloads behave very linearly for reasonable 
airfoil motions. Reference [51] included the static twisting of the 
airfoil due to the steady pitching moment into the flutter problem 
and demonstrated a marked effect upon flutter boundaries, par- 
ticularly for the supercritical MBB-A3 airfoil. Reference [28] 
documented the utility of s-plane Pade curve fits of transonic air- 
loads (which rely on the concepts of linearity and superposition) 
for aeroelastic analysis. Viscous effects are shown to generally 
result in larger values of flutter speed since transonic effects are 
alleviated by the boundary layer. Zimmerman [170] describes 
the application of several unsteady aerodynamic theories to the 
prediction of transonic flutter instabilities for a 2D airfoil pro- 
file. Flutter predictions are made using TSP theory, linearized 
TSP, TSP with quasi-unsteady boundary-layer and doublet lat- 
tice. Both frequency domain and time-marching methods are 
used in these flutter predictions. 

Wu et al. [164] describe the use of a compressible-flow, 
FNS code for the prediction of unsteady airloads and flutter 
of 2D airfoil sections. Unsteady flow results are shown for 
harmonically oscillating NACA 0006 and NACA 0012 airfoils. 
Time marching aeroelastic solutions are given for the NACA 
0012 airfoil showing stall flutter for several different Mach 
numbers and Reynolds numbers. 

At the NLR, a number of unsteady aerodynamic computa- 
tional procedures have been developed for use in transonic flut- 
ter prediction (Hounjet [74] and Steiginga and Houwink [149]). 
Zwann [168] describes an investigation into the accuracy of 
several of these methods for predicting flutter of a transport- 
type wing model. The methods examined include doublet lat- 
tice with quasi-steady corrections, the Q3D quasi 3D method 
and the FTRAN3 hybrid field-panel/finite-difference procedure. 
Corrected doublet lattice results are shown to display a flutter 
dip, probably due to the use of available experimental sectional 
airloads data. The flutter dip computed is, however, uncon- 
servative. The correlations of flutter speed versus Mach num- 
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ber reported for the Q3D code are slightly conservative, but in 
good agreement with test results over a range of Mach numbers, 
and up to a point where the flow begins to separate. Several 
FTRAN3 predicted flutter points compared closely to the Q3D 
results. Hounjet and Meijer [75] give additional applications of 
the time linearized F I RAN 3 code to flutter calculations for a 
fighter configuration with stores. 

Other time-linearized finite-difference methods for solving 
the unsteady transonic flow about harmonically oscillating wings 
are described by Weatherill and Ehlers [161] and Shieh, Schoen, 
and Fung [143]. In Reference [161] the transonic small distur- 
bance equation is split into a steady and an unsteady flow com- 
ponent, and the unsteady portion is solved using an out-of-core 
LU decomposition procedure (OPTRAN3). Sample flutter anal- 
ysis results are presented for a flexible rectangular-wing. Refer- 
ence [143] describes an implicit finite-difference procedure for 
solving the 3D full-potential equation for small unsteady per- 
turbations of the body surface. The potential function is split 
into steady and unsteady components. The steady component 
can be obtained from a number of sources, such as higher order 
fluid dynamic equations (e.g. Euler, Navier-Stokes) or experi- 
mental data. Corrections are applied for nonisentropic effects, 
and the body motion is enforced using surface transpiration on 
a static grid conforming to the mean position of the solid sur- 
faces (ZUNAS code). Computed results are given for unsteady 
oscillations of a NACA 0012 airfoil and for flutter boundaries 
of a 3D transport-type wing. 

Comparisons of calculated and experimental flutter bound- 
aries for wings have been given by Guruswamy and Goorjian 
[66], Isogai [87], Isogai and Suetsugu [88] and Myers et al. 
[113]. Isogai [85] studied the supercritical wing of Farmer et 
al. [57] using a nonconservative full potential code and an inter- 
acted boundary layer model. The trend of the transonic flutter 
dip is very nicely predicted although the dip occurs about 0.08 
low in Mach number. The premature flutter dip and the subse- 
quent premature rise of the calculated boundary at higher speeds 
is of concern. Myers [113] also shows such a premature rise 
in the boundary. Finally, Isogai [87] shows flutter comparisons 
for a different supercritical transport wing which agree nicely 
with the experimental flutter dip. 

Cunningham et al. [44] describe TSD code flutter calcula- 
tions for the AGARD Standard Aeroelastic flutter model con- 
figuration [166]. This is a 45° sweptback wing with a taper 
ratio of 0.66. Experimental and computed flutter boundaries are 
given for Mach numbers from 0.338 to 1.14. For this 4 per- 
cent thick wing, transonic effects are delayed to high subsonic 
Mach numbers and linear theory results from both CAP-TSD 
and a kernel function program are in very good agreement with 
experiment up to Moo = 0.98. Nonlinear CAP-TSD subsonic 
flutter calculations agree better with experiment than the linear 
theory, particularly for the change in slope of the flutter bound- 
ary near M 0 0 = 0.95. Robinson et al. [130] have also performed 
flutter calculations for this wing using an Euler code. The mov- 
ing mesh scheme described above was modified for use with 
a structured grid algorithm and the flutter boundary prediction 
was in good agreement with the TSD code and experiment. 

Borland and Nagaraja [37] describe extensions to the 
XTRAN3S code to permit aerodynamic modeling of wing- 
pylon-stores and wing-fuselage configurations as well as the 
addition of supersonic analysis capability, non-reflecting bound- 


ary conditions and a Fourier analysis option for unsteady surface 
pressures. A number of applications of the code are presented 
for steady- and unsteady-flow aerodynamics and static/dynamic 
aeroelastic solutions. The configurations studied included an A- 
6 fighter configuration with stores, and a transport wing/fuselage 
configuration with winglet. 

Mulak, Meurzec and Angelini [112] describe a finite- 
difference procedure to solve the 3D unsteady transonic small 
disturbance equation. The method solves a discretized form of 
the equation using an ADI splitting technique. Unsteady gener- 
alized airforces are computed for a number of normal modes for 
a swept, high-aspect-ratio transport-type wing. Some compar- 
isons are made for steady and unsteady surface pressures and 
test data at a transonic Mach number. Flutter calculations are 
performed using different wing twist distributions correspond- 
ing to jig and flight condition shapes. A comparison is also 
presented for flutter calculations using both the doublet lattice 
and the TSD method. 

Wong and Lee [165] describe a numerical procedure to com- 
pute unsteady transonic flow about 3D isolated-wing configura- 
tions (UST3D code). The procedure solves the transonic small 
disturbance equation, split into a nonlinear steady component, 
and a linear unsteady component. The steady part of the equa- 
tion is solved with an implicit Newton-iteration scheme, used in 
conjunction with a preconditioned gradient type procedure. The 
unsteady part is solved with a semi-implicit technique which is 
explicit in the stream direction and implicit in the cross-flow 
planes. A nonreflecting far-field boundary condition is applied 
at the outer boundaries of the Cartesian mesh. Computed results 
are given for a fighter wing. Real and imaginary components 
of lifting pressure are given for a Mach number of 0.9 and k = 
0.45. Aeroelastic results are compared with doublet lattice com- 
putations for unsteady pressures and a sample flutter calculation 
is presented. 

Five additional applications of TSD codes to wing flutter 
calculations are described in Ref. [13], Bennett and Batina [26], 
Gibbons et al. [62], Silva and Bennett [144], and Guruswamy 
and Gooijian [69]. Reference [13] describes an implicit finite- 
difference procedure (ATRANS code) for solving the 3D tran- 
sonic small disturbance equations on a Cartesian mesh. Aeroe- 
lastic calculations are made for the Falcon 900 wing at tran- 
sonic Mach numbers, both with and without the incorporation 
of static aeroelastic deformations. Comparisons are also made 
with a doublet lattice method at a subsonic Mach number of 0.8. 
Bennett and Batina [26] present flutter calculations for a three 
percent thick clipped tip delta wing with a leading edge sweep 
angle of 50.5 degrees. The Mach range covered was 0.6-0.9. 
At the lower Mach numbers the results were in excellent agree- 
ment with linear theory, while there was a 6 percent reduction 
in flutter speed from linear theory at Moo = 0.907, bringing the 
results into better agreement with experiment. Gibbons studied 
a 70 degree swept delta wing for Mach numbers ranging from 
0.6 to 3.0. The calculated flutter speeds are in very good agree- 
ment with experiment at transonic speeds. At supersonic speeds, 
the trend of the flutter boundary with Mach number agrees well 
with experiment. Silva and Bennett [144] show transonic flutter 
boundary predictions for a complex wind tunnel model. The 
blended wing-fuselage and tip ballast store were modeled. Sig- 
nificant improvement in predicted flutter points over those cal- 
culated with linear theory are shown for high transonic Mach 
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number experimental points. Finally, Guruswamy and Gootjian 
[69] give flutter boundary calculations for a rectangular wing 
which agree well with experiment including a low supersonic 
condition. 

Several investigators have studied airfoil and wing limit- 
cycle-oscillations (LCO) using both CFD-based and empirically- 
based computational methods. Bendiksen and Kousen [25] and 
Kousen and Bendiksen [94, 95] studied nonlinear aeroelastic 
dynamic response of an airfoil section using a CFD method 
based on the Euler equations. The Euler procedure is a finite- 
volume scheme solved using an explicit five-step Runge-Kutta 
algorithm. A two-degree-of-ffeedom structural dynamics model 
was solved iteratively along with the fluid dynamic equations 
to predict the dynamic response of the airfoil section. Aeroe- 
lastic stability was studied for NACA 0006 and NACA 64 AO 10 
airfoil sections over a range of Mach numbers. Limit cycle 
oscillations were found to occur at transonic speeds as a bifur- 
cation phenomenon. A divergence/flutter interaction [94] was 
identified for which the airfoil dynamically approaches a static 
equilibrium offset at non-zero angles-of-attack. 

Kousen and Bendiksen [95] also examined the effects of 
torsional spring free-play on the flutter response of 2D airfoil 
sections restrained by a two-degree-of-freedom elastic system 
(pitch and plunge). For the 2D case, the incorporation of free- 
play in the torsional spring lowers the stability boundary for 
conventional flutter. A bifurcation phenomenon is shown to ex- 
ist above the maximum reduced velocity for flutter. Above these 
values, a limit-cycle behavior is exhibited by the airfoil/spring 
system. This 2D analysis method was applied to the flutter be- 
havior of a well known typical section model of sweptback wing 
bending and torsion modes. It is shown that the dual-mode in- 
stability calculated for this model corresponds to a region, or 
pocket, of limit-cycle behavior. Finally, the effects of reduced 
airfoil thickness were examined and shown to be destabilizing, 
leading to limit-cycle behavior at lower Mach numbers than for 
a thicker airfoil section. 

Meijer and Zwann [110] describe a semi-empirical compu- 
tational method for the prediction of LCO occurring on fighter- 
type aircraft (see figure 5). Steady-state experimental surface 
pressure data are used to compute generalized airforces required 
to solve the structural dynamic equations of motion for a ve- 
hicle. The surface pressures are tabulated versus Mach num- 
ber and steady angle-of-attack and utilized during the transient 
response simulation in a quasi-steady fashion wherein the in- 
stantaneous, induced angle of attack is used to reference the 
pressures. Both rigid and elastic mode shapes are used in the 
computations. Sample calculations of dynamic structural re- 
sponses were made for two vehicle configurations, an isolated 
wing and a wing-with-missile configuration. Limit-cycle oscil- 
lations were found only for the clean-wing configuration, where 
the computed frequency and mean angle-of-attack data com- 
pared well with flight test results. For the flight conditions were 
LCO was observed, single-degree-of-freedom predictions using 
single mode shapes failed to show indications of LCO behavior. 
The authors state that LCO may not be associated with single 
degree-of-freedom systems, but rather with multiple-degree-of- 
freedom systems where modal coupling takes place, in a manner 
similar to the classical flutter mechanism. 

The suppression of flutter instabilities by active controls has 
been demonstrated using CFD methods by Ominsky and Ide 


[119] and Guruswamy [68]. Ominsky and Ide use a Full Po- 
tential equation CFD method to evaluate aeroservoelastic con- 
trol laws for flutter suppression. A test case for a swept wing 
is created by first determining a flutter condition to be con- 
trolled. A single trailing edge flap is deflected with time using 
a simple open loop multi-parameter control law. The param- 
eters of flap amplitude of oscillation frequency and phase lag 
are studied to determine their effects on flutter response and 
a second-order control law is used in a closed loop manner to 
control the identified flutter instability. Guruswamy summarizes 
efforts at NASA-Ames to develop CFD methods for aeroelastic 
applications. Applications of the ATRAN3S code (TSD) and 
the TNS code (Euler/Navier-Stokes) are given. These include 
unsteady aerodynamic results, flutter analyses and flutter sup- 
pression demonstrations. Some comparisons with experimental 
data are given. 

Schuster et al. [139] have developed a FNS code using 
zonal grid generation techniques to enable aeroelastic analysis of 
complete vehicles. Static aeroelastic calculations are presented 
for a fighter aircraft at a high wing-loading transonic condition. 
Calculated twisting of the aeroelastically tailored wing/fuselage 
configuration compares favorably with experiment. 

Obayashi and Guruswamy [118] describe applications of the 
ENSAERO code to compute unsteady aerodynamics of a delta 
wing in oscillatory and ramp motions. The ENSAERO code 
is an Euler/Navier-Stokes method which includes a structural 
dynamics model to account for wing flexibility. The method 
solves the thin-layer Navier-Stokes equations and incorporates 
an algebraic turbulence model. Comparisons of theory and ex- 
perimental data are presented for a number of steady flow and 
unsteady flow (pitching-oscillations) conditions. Computational 
results are also presented for a delta-wing undergoing an un- 
steady ramp motion (pitch). The computed results indicate a 
vortex breakdown condition occurring for a ramp motion which 
reaches a maximum amplitude of 10 degrees angle-of-attack. 

Vinh et al. [159] present an interesting augmentation to the 
aeroelastic response capability of such codes. In addition to the 
structural dynamic equations of motion normally used for aeroe- 
lastic analysis, two additional degrees of freedom modeling the 
rigid aircraft short period mode were implemented along with 
an automatic trimming capability. This allows the interaction 
of the elastic modes upon vehicle stability and control to be 
calculated in a straightforward manner. 

7.4 Viscous-Inviscid Int eraction Applications 

The transonic aeroelastic stability illustrations shown in fig- 
ures 3-5 all involve strong viscous effects which inviscid analy- 
ses cannot predict. In order to achieve desired improved accu- 
racy in predictive methods for such cases, an accounting of the 
effect of the viscous boundary layer is mandatory. Rows which 
must be treated include those with thickening boundary layers, 
incipient separation, shock-induced and trailing-edge separation, 
and alternately separating and reattaching flows. The computa- 
tional resource estimates given in Table II indicate the severe 
cost incurred by going to higher level CFD flow models and 
lead to the desirability of achieving the maximum possible vis- 
cous flow analysis capability with lower level CFD methods 
coupled with interacted viscous boundary layer models. Due to 
the importance of this issue for aeroelastic analysis, this sec- 
tion focuses on such applications. Interacted boundary layer 
methods have been developed for unsteady applications using 
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direct solution for attached flows and indirect solution for sep- 
arated flows. Calculations using TSD, full potential, and Euler 
equation codes have been extensively reported. 

Examples of viscous flow effects for attached 2-D flows 
(direct solution method) are given by Guruswamy and Goorjian 
[65], Howlett [80, 81], and Houwink [77, 78]. All of these 
references use an integral boundary layer model coupled to a 
TSD inviscid method and calculations for the AGARD Stan- 
dard Configurations are given (NACA 64 AO 10, NACA 64A006, 
MBB-A3, NLR 7301, and NACA 0012 airfoils). Pirzadeh and 
Whitfield [122] report a 3-D direct viscous solution coupled 
to the Euler equations and give results for the NACA 64 A0 10 
airfoil and the ONER A M6 wing. 

Unsteady results using the quasi-simultaneous method cou- 
pled to TSD solvers are reported by Houwink and Veldman [79] 
for 2-D cases and by Henke et al. [71] for 3-D cases. Houwink 
gives separated flow results for an oscillating supercritical airfoil 
and for an airfoil with a deflected spoiler. Henke gives details 
of a similar viscous coupling method implemented in a stripwise 
fashion in a 3-D TSD code and including aeroelastic compu- 
tations. Comparisons with measured unsteady pressure from 
a modem transport wing model with a supercritical airfoil are 
shown and flutter calculations using the interacted viscous-TSD 
code are compared with doublet lattice calculations. 

The most elaborate viscous interaction calculations of un- 
steady separated 2-D flows have been published by LeBalleur 
and his coworkers [98, 99]. An unsteady deflect integral bound- 
ary layer model is coupled to a TSD solver using the alternating- 
direction implicit solution method which enables an incorpora- 
tion of the strong interaction of the viscous influence during the 
final z-sweep. The semi-implicit solution procedure involves a 
relaxation solution of a viscous influence function to force con- 
vergence of the viscous and inviscid solutions at each time step. 
Reference [99] gives unsteady shock-induced separation results 
for an oscillating NACA 64A010 airfoil, self-induced shock os- 
cillations for a circular arc airfoil, and a supercritical airfoil with 
an oscillating spoiler. Giroudrouux-Lavigne and LeBalleur [63] 
further explore the self-excited shock oscillations, giving exam- 
ples for the RA16 SCI supercritical airfoil and the NACA 0012 
airfoil. The shock oscillation conditions are shown to coincide 
with experimentally observed buffet onset boundaries. It is in- 
teresting to note similar calculations of buffet onset boundaries 
reported by Hirose and Miwa [72]. Experimental buffet bound- 
aries for the NACA 0012 and KORN 75-06-12 airfoils are 
compared with calculations from a 2-D TNS code. As in the 
above viscous-inviscid interaction results, self-excited oscillat- 
ing shocks are observed for conditions in good agreement with 
the experimental buffet boundaries. 

In summary, viscous interaction methods are capable of 
treating important transonic effects when coupled with lower 
level CFD methods. Impressive results have been published 
for 2-D flows. Similar capability for 3-D flows involving 
shock-induced separating and reattaching flows remains a most 
important research topic. 

7.5 Complete Aircraft Configurations 


ered are for steady conditions; thus they can serve as a baseline 
for estimating resources which would be required for unsteady 
analyses. It is not surprising that the most detailed aircraft ge- 
ometry modeling has been applied for cases at low angles of 
attack, near design conditions, where the flow is attached or 
mildly separated. These conditions are most important in de- 
sign and the codes can be expected to perform at their best due 
to good flow quality (steady, attached, thin shear-layers, etc.). 
Thus, indications of the accuracy of results here will help to as- 
sess the readiness of the codes for the more demanding transonic 
flutter and low-speed, high-angle buffet conditions. Table III 
summarizes a number of CFD applications for such cases. All 
of these studies used the TLNS equations and all implemented 
simple turbulence models. All used convergence acceleration 
devices (local time-stepping, multigrid, etc.) which yield accu- 
rate results only for converged, steady flows. Except for Flores 
and Chaderjian [59], all made use of assumed symmetry in the 
flow to reduce by half the grid size (only one-half aircraft model 
analyzed). Finally, the codes used in these studies were gener- 
ally second-order accurate in space and first-order accurate in 
time. Fujii and Obayashi [61] and Huband et al. [83] give 
results for complete aircraft models at transonic speeds while 
Thomas et al. [151] studies a detailed forebody/strake model 
at low speed and high angle. The latter is included due to its 
complex surface modeling. Vadyak and Schuster [156] give re- 
sults for a sharp-edged generic fighter wing-body configuration 
at high angles. Kwon and Sankar [97] and Rizk and Gee [128] 

give results of buffet flowfield calculations about two high per- 
formance fighter configurations at high angles. These last three 
cases of high angle vortex-dominated flows will be discussed 
in the next section. 

Fujii and Obayashi [61] modeled the W-18 transport con- 
figuration as a wing-fuselage and made calculations for three 
angles-of-attack. The overall surface pressures compared fairly 
well with experiment except in the outboard wing region where 
elastic deformations were not accounted for. Their numerical 
algorithm allowed the calculations with 700K grid points to be 
obtained in 5-6 cpu hours. Flores and Chaderjian’s [59] study 
of the F-16A aircraft is one of the most ambitious applications 
to date. They modeled the wing, fuselage, tail, inlet, inlet- 
diverter, and the exhaust nozzle using 27 grid zones for the 
half-airplane. The flow-through inlet was modeled including 
power effects. Comparisons with experimental pressures indi- 
cated that the wing leading-edge expansion was not adequately 
resolved and the wing shock location was off by 12% chord. 
Doubling the grid size to one million points allowed the cal- 
culation of the full-span aircraft at five degrees sideslip angle. 
This also doubled the runtime to 50 hours. With regard to buf- 
fet calculations, it is interesting to note that the vertical tail tip 
vortex for this condition dissipated within one tip chordlength 
due to numerical dissipation (due of grid stretching downstream 
of the tip). Huband et al. [83] studied the same F-16A (the inlet 
is faired over) for a low supersonic Mach number. Their fine 
grid solution (1,24 IK grid points) occupied 59 million words of 
memory and required 40 hours of runtime. They obtained fa- 
vorable agreement with experimental surface pressures but the 
wing leading-edge suction peaks were not correctly predicted 
due to lack of numerical resolution. 


This section will survey the most ambitious CFD modeling 
of complete aircraft in order to indicate the levels of details and 
the level of effort such studies involve. All of the cases cov- 
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TABLE m — CFD APPLICATIONS TO COMPLEX AIRCRAFT GEOMETRIES 


Ref. 

Configuration 

M 

a, deg. 

Re, 10 4 

N« 

Njp 

r, fiscc 

Tq, u » hrs j 

i 

Fujii 41 

. 

W-18 

0.82 

2.4,6 

1.67 

300G+ 

700K 

9.5 

5-6 | 

Flores 59 j F-16A w/inlet 

0.9 

6 

4.5 

5000 

500K 

36+ 

25 i 

i 

Huband 83 F-16A w/o inlet 

L 2 

6 

12.75 

40,000+ 

U41K 

30+ 

40 

Thomas 131 

F/A-18 forebody 
& LEX 

0.3 

30 

0.74 WT 
10.0 FLT 


300K 


: 

Vadjak 154 | generic wing/body 

i 

0.3 

10,21*. 

23,25 

1.25 j U00 

i 

101K 


i 

Kwon 97 S F-15 wing/body/inlet 

i 

0.15 

20* 

1 

224K 



^^128 | F-18 aircraft 

0.243 

30* 

11.0 

4,100 



55 


♦Unsteady vortex flows 


TABLE IV — CFD CALCULATIONS FOR HIGHLY-SWEPT WINGS AT HIGH ANGLES OF ATTACK 


Ref. 

Configuration 

A, deg. 

q, deg. 

Re, 10 6 

Moo 

Ngp,K 

N st 

T,/i sec 

Tcp U ,hrs 

Hitzel 73 

cropped delta, 
s,e 

65 

10,20,24,28* 

10-28,24*, 

28*10 

1.25 

0.4 

0.85 

1.20 

270, 540 




Hsu 82 

delta, r 

60 

12(20) 

1.3 

0 

712 

700(1000) 

52 

7 (10)+ 

Hsu 82 

double-delta, r 

80-60 

12(20) 

0.01, 
1.3, 100 

0 

955 

700(1000) 

52 

10 (14)+ 

Hsu 82 

cropped, 
double -delta, r, 

80-40 

12(20) 

1.3 

0 

1021 

700(1000) 

52 

10 (14)+ 

Hsu 82 

double-delta, r 

80-60 

6,12,15,20, 

25,30,35*, 

40* 

1.3 

0 

859 

500(1000) 

52 

6 (12)+ 

Fujii 60 

double-delta, r 

80-60 

6,12,(30*,35*) 

1.3 

0.3 

850 

1000(5000) 

20 

5 (25)+ 

Thomas 150 

delta, s 

75 

0-40, 

20.5,40* 

0.95 

0.3 

545 

400 

40 

2.5 

Kris: 96 

delta, s 

75 

20.5 

0.5 

0.3 

350(1000) 




Agarwal 1 

delta, s 

70 

(20, 25, 30, 
35, 40)* 

1.0 

0.3 


1-2000 

. 


3 


e-Euler code 


s-sharp leading-edge 


r-round leading-edge 


‘-vortex burst/unsteady calculations 


•(-estimated 
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8. High-Angle Vortex-Dominated Flow Applications 

Computational studies are beginning to delineate the re- 
quired levels of effort to produce accurate computations for sep- 
arated, vortex-dominated flows. In this section attention will be 
given to a brief assessment of the state-of-the-art in computing 
vortical flows using CFD methods, followed by a discussion of 
recent research efforts in two specialized categories, wing rock 
and vortex buffet. 

8.1 Status of CFD Prediction Methods for Vortical Flows 

Edwards [56] provides a detailed study of current research 
and applications into the use of CFD methods for the predic- 
tion of separated, vortical flows over airfoils, wings and bodies. 
Table IV lists a selection of these applications relevant to com- 
putational aeroelasticity at high angles. All except Hitzel [73] 
represent applications of the TLNS equation, while Hitzel uti- 
lizes the Euler equations. The major conclusions of this study 
are repeated here for completeness. 

Calculations of stable vortex -dominated flows about simpli- 
fied highly swept wing geometries at low speeds are available. 
Early Euler equation results gave encouragement to their use for 
such flows, but recent detailed studies point to issues regarding 
their accuracy and convergence for realistic geometries. The 
thin-layer Navier-Stokes equations, both laminar and turbulent 
are being used. The results capture key features of the flowfield 
such as primary and secondary vortices and surface pressure 
details. Indications that vortex breakdown is being simulated 
have been published. However, all cases surveyed indicated the 
need for further grid refinement in order to achieve quantitative 
agreement with experiment. Currently, grid sizes of 250,000 to 
10 6 grid points are being used for half airplane modeling lead- 
ing to CPU runtimes of 2-25 hours. These grids are typical of 
those developed for attached shear flows with the highest grid 
density near the body. This leads to inadequate density in off- 
the-body regions where concentrated vortices are located and 
mesh enrichment methods are being developed to address this 
problem. A fundamental problem is the lack of any turbulence 
model designed for dynamic free shear layers. 

Agrawal et al. [1] provide an interesting numerical study of 
this issue of off-the-surface turbulence modeling. Calculations 
for a sharp-edged delta wing from an Euler code, a laminar 
TNS code and a turbulent TNS code show vortex bursting at 
locations in general agreement with experiment. It is noted that 
the solutions are based on meshes that are not usually considered 
fine enough for resolving flowfields in primary vortex regions. 

8.2 Rolling Oscillations and Wing Rock 

Several studies have been published of unsteady airloads on 
rolling delta wings at high incidences. The dynamics of the 
leading-edge vortices which form at such conditions are one 
of the driving forces which can lead to wing-rock; a limited 
amplitude rolling oscillation which limits aircraft performance. 
Computations have been reported for low speed flows and for su- 
personic flows simulated via a conical flow assumption. Nayfeh 
et al. [116] use an unsteady vortex-lattice method to predict 
low-speed experimental wing-rock conditions. 


Kandil and his coauthors have studied unsteady vortex- 
dominated flow about delta wings using the conical flow as- 
sumption of supersonic flow. Reference [90] gives compress- 
ible Navier-Stokes results for rolling oscillations while Ref. [91] 
gives results for oscillating leading-edge flaps. Lee and Batina 
[101], again making use of the conical flow assumption, couple 
an unstructured grid Euler code with a wing rolling equation 
of motion to make calculations of supersonic wing-rocking mo- 
tions. In a later work, Lee and Batina [100] demonstrated that 
this supersonic wing-rock phenomenon can be controlled with 
the use of oscillating, leading-edge flaps, activated by a simple 
control-law. More recently, Kandil and Salman [92] used an 
Euler conical flow solution procedure to study the control of 
wing-rock using leading edge flaps, together with a nonlinear 
control law. 

8.3 Buffet Flowfield Calculations 

Current and future military fighters are called upon to per- 
form high angle maneuvers at elevated loading conditions and 
are experiencing structural fatigue problems due to dynamic buf- 
fet loads on aft components of the empennage, in particular the 
vertical tail(s). This has led to a desire for improved prediction 
methods for such flows, which emerging CFD methods might 
be expected to fulfill. Key flow modeling issues involve: i) un- 
steady flow separation for moderately swept, rounded leading- 
edge wings, ii) grid densities necessary for accurate calculation 
of vorticity convected over significant distances, iii) dynamic 
turbulence modeling for free shear layers, iv) detailed modeling 
necessary for accurate calculations of vortex instabilities, and v) 
vortex flows about wings at high subsonic speeds. 

Three cases of the computation of buffet flowfields are given 
as the last three entries in Table III. All are for quite complex 
configuration models of high performance fighter aircraft and it 
has been noted above that this type of computation is the most 
demanding for aeroelastic. applications. Vadyak and Schuster 
[156] made calculations for a generic fighter configuration con- 
sisting of a sharp-edged strake/wing/fuselage. Low speed wind 
tunnel LDV flowfield data was available and crossflow velocity 
comparisons for a = 21 deg. appear to be good. A bubble-type 
reverse flow region, indicating breakdown, is observed at ap- 
proximately x/L = 0.9. These comparisons for this high angle, 
vortex flow case are very noteworthy, particularly due to the 
complex geometry tested. However, the computational results, 
obtained with a grid of 101,000 points, cannot be regarded as 
converged. 

Kwon and Sankar [97] give calculations for a half- span 
wing/body/inlet model of the F- 15 aircraft. The FNS equations 
are solved using a hybrid time differencing scheme suitable 
for implementation on virtual memory machines. On a grid 
consisting of 224,000 points, calculations of unsteady, buffet- 
like flowfields are given for M oo =0.15 and a = 20°. The 
unsteadiness is shown to consist primarily of a low frequency 
of about 5 Hz. (full scale aircraft) and a high frequency 
ranging from 29-34 Hz. The higher frequency content compares 
favorably with low speed wind tunnel tests for the complete 
aircraft (Triplett [154]) even with the relatively coarse grid used. 

The final entry in Table III, that of Rizk and Gee [128], 
gives results for a complex configuration model of the half- 
span F-18 aircraft. An overset zonal grid scheme consisting 
of ten grids was used to model the forebody, fuselage, LEX, 



17 


faired-over inlet, wing, deflected leading-edge flaps, vertical and 
horizontal tails and an idealized boundary layer diverter vent. 
TLNS results are given for A/**, = 0.243 and a = 30.3°. The 
solution required about 8 million words of memory and took 
about 55 hours of CPU time for 5,100 steps. Bursting of the 
LEX vortex in the vicinity observed in flight tests is shown. 
Unsteady loads on the vertical tail show a dominant frequency 
of 15-20 Hz. Other higher frequencies are evident but are not 
well resolved. Triplett [154] gives test results indicating that the 
principal modal buffeting response frequencies of this aircraft 
are 15 Hz. and 45 Hz. Earlier calculations [127] indicate 
that finer meshes (e.g. about 2 million grid points) will be 
needed to accurately predict such features as the location of 
vortex bursting. 

The buffeting response of launch vehicles is a critical issue 
in their performance. Up to the present, wind tunnel studies 
have largely been relied upon for the prediction of such airloads. 
Azevedo [14] presents novel TLNS results of buffeting response 
calculations for a hammerhead payload configuration. 

In summary, ambitious applications of CFD to complicated 
aircraft geometries are being performed. Available studies are 
very encouraging in the overall agreement with experiment. Im- 
portant areas are also being highlighted where additional grid 
resolution is needed to achieve local accuracy in such impor- 
tant features as suction peaks and shock locations. Computer 
runtimes for these cases are in the range of 5-55 hours (half- 
airplane). It is probable that accurate buffet calculations will 
require capability such as this forebody and wing flows in or- 
der to generate accurate “starting” conditions for the convecting 
vortex flows. It is also probable that similar capability will be 
required to calculate the buffeting response of the tail to the on- 
coming buffet flow. Indeed, this is likely to be a more difficult 
problem than the calculation of the “starting” wing flow due to 
the turbulent nature of the local flow at the tail. 

9. Discussion and Future Trends 

The proceeding sections have surveyed a segment of the 
large efforts that have been spent on computational unsteady 
aerodynamics suitable for aeroelastic applications of fixed wing 
aircraft. These efforts can be divided into five general cate- 
gories: i) development and demonstration of unsteady aerody- 
namic methods, ii) rigorous calibration and validation of these 
unsteady methods, iii) development and demonstration of aeroe- 
lastic analysis procedures, iv) rigorous calibration and validation 
of these aeroelastic procedures, and v) application of the result- 
ing aeroelastic methods in the design or modification of actual or 
proposed flight vehicle configurations. Although nearly a decade 
and a half has passed since the earliest pioneering demonstra- 
tions of CFD applied to aeroelastic analyses, considerable efforts 
are still needed to reach a state where the aircraft industry, as a 
whole, accepts and routinely utilizes CFD for category v. activ- 
ities. By far, the largest number of references cited herein fall 
within category i., and to a lesser extent, category iii. Also, to 
a large extent, the authorship of these cited references seems to 
indicate that a large segment of this research is currently occur- 
ring in universities and government sponsored research labora- 
tories. Although there are some notable exceptions to this trend, 
many more industrial applications of aeroelastic CFD methods 
are needed and desirable. The following comments are offered 
regarding the current status and future trends of computational 
methods to meet the challenges listed in Table I and illustrated 


in figures 3-7: 

1. Stability and control estimates, including aeroelastic effects, 
are being treated with linear methods which are well suited 
for design conditions well removed from flow separation 
boundaries where viscous effects are important. Initial 
demonstrations of the unified treatment of complete vehicle 
dynamics (rigid body and aeroelastic) using CFD methods 
are available. 

2. Gust response analyses are currently conducted using linear 
aerodynamic methods. While CFD codes capable of treating 
complete configurations (necessary for gust response anal- 
ysis) are available, there have been no reported attempts of 
such calculations other than isolated vortex-airfoil interac- 
tions. 

3. A number of studies of flutter boundary calculations with 
CFD codes have been cited. Many more cases are needed 
in order to establish the required level of flow modeling de- 
tail which will provide accurate and useful results. It is not 
yet established for what conditions nonlinear flow modeling 
results differ sufficiently from linear results to justify their 
use. Also, the tradeoffs between the computational expense 
of the various levels of flow modeling (e.g. TSD-viscous 
vs. TLNS) and improved accuracy of results have not yet 
been established. A few cases have shown significant im- 
provement over linear theory for small, but critical, ranges 
of transonic Mach numbers. 

4. Prediction of these critical minimum transonic flutter speeds 
and nonclassical aeroelastic response phenomena, such as 
control surface buzz and aeroelastic limit cycle oscillations, 
certainly will require, at the very least, reliable, robust 
interacted boundary layer models capable of handling some 
amount of flow separation and reattachment. This capability 
is not yet mature for wings or more complete configuration 
modeling. This is also the case for the important design 
issues relating to aerodynamic control surfaces: control 
effectiveness and control hinge moments. 

5. Cases of tail buffeting and structural fatigue being encoun- 
tered by operational aircraft have focused attention upon this 
area. The achievements in the ability to compute vortex- 
dominated flows are truly impressive and initial calculations 
of buffet-like flows appear to contain elements (frequency 
content, etc.) of realism. However, studies of the accuracy 
and convergence of such calculations in terms of grid den- 
sity, sensitivity to initial conditions, repeatability, etc. are 
necessary. Two interesting questions regarding such flow- 
field simulations are: 

a. if the flowfield simulations are repeatable, can the flow 
modeling contain adequate fluid dynamic modeling to 
ensure accurate buffet flows calculations? 

b. if the flowfield simulations are very sensitive to initial 
conditions or are not repeatable, what measures are to 
be used in establishing the required number of cases 
and computational record lengths to ensure solution 
convergence? 
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Regarding future trends in computational aeroelasticity the fol- 
lowing comments are offered: 

1. Robust interactive boundary layer methods coupled with all 
levels of inviscid codes will mature and be applied to aeroe- 
lastic stability and response calculations for high-speed, low 
angle conditions. The resulting improved accuracy in key 
areas such as transonic flutter prediction will be of great 
value. 

2. Studies will delineate conditions best suited for the use of 
CFD methods. A hierarchy of the problem areas, keyed 
to the complexity of the flow problem being studied (e.g. 
Types I, II, or III flows), will be established and will be 
addressed by a suite of tools drawn from the alternative 
levels of CFD flow modeling. 

3. Unstructured and hybrid grid methods will mature, leading 
to computations of quite detailed configuration complex- 
ity. The increased expense of unstructured grid methods 
(due to more complex coding) will be offset by the antic- 
ipated payoffs of the approach: fewer grid points required 
for comparable accuracy along with the complex geometry 
modeling capability. Continued improvements in algorithm 
efficiency will be seen for all gridding approaches; block- 

structured, zonal and hybrid grid schemes will also be used 
in the search for robust, economical and user-friendly ge- 
ometry modeling capability. 

4. Improvements in conventional computer architectures will 
not achieve the orders of magnitude speedup needed to make 
the most complex aeroelastic computations, such as tail buf- 
feting, routine. Massively parallel processing architectures 
may eventually provide relief, but not in the near future. 

5. The need for experimental data bases specialized for the val- 
idation of computational methods for interesting aeroelastic 
flows is being addressed by NASA’s Benchmark Models 
Program and the European Aeroelastic Models Program. 
Tests directed at unsteady off- the -surface flow measure- 
ments will continue to press the limits of instrumentation 
technology. 

6. A final comment relates to differing uses of computational 
methods in the fields of aircraft performance analysis and 
aeroelastic analysis. The former application may be broadly 
typified as dealing with well-behaved, stable flows exempli- 
fied by performance and design studies where steady CFD 
methods are currently being utilized heavily. The latter 
application typically must deal with flows on the edge of 
flow breakdown (i.e. structural dynamic instabilities in the 
presence of dynamic flows and flows transitioning between 
states) where safety-of-flight is of paramount concern. At- 
tention by the user-communities to the appropriately differ- 
ent standards to be applied to the validation or calibration 
of CFD methods in these two areas is required. 

In conclusion, it is anticipated that unsteady CFD methods 

will emerge to be routinely used for aeroelasucally challenging 

problems. 
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